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Abstract 

This work addresses the topic of knotted stable structures in excitable media. These structures 
appear in a wide variety of situations, such as cardiac fibrillation, chemical reactions, etc. Entangled 
curves have been found in numerical computations of the equations that describe excitable media. 
They present an unusual stability. An explanation for this behaviour has been an open question. 
In the present work we introduce for the first time the meaning of the helicity in an excitable 
media as a new tool to study the stability of these systems. The helicity is related to the total 
entanglement of the system. We have studied how the helicity is conserved or lost through the 
walls of the medium and shown that these behaviours are dominated by the boundary conditions, 
so the distortion of these conditions could lead to the disappearance of the structures. 

PACS numbers: 82.40. Ck, 47.32. C, 47.32.cd 
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There has been experimental reports on stable structures (rotors) in cardiac muscle that 
appear in cardiac fibrillation which may support the ideas of Arthur T. Winfree {2] on 
the relation between electrical wave propagation in the heart and some stable structures 
appeared in excitable media. One of the most interesting aspects of excitable media is 
that they support vortices, called spiral or scroll waves. Spiral waves are characterized 
by the fact that they rotate around a topological point defect called phase singularity or 
rotor. Since spiral waves are the cross sections of a scroll wave, the phase singularity can be 
interpreted as the cross section of a filament in three dimensions. These curves can be linked 
or constitute knots. The important fact of the phase singularities is that they are the focus 
or organizing centres of the excitation of the complete medium, and that their stability has 
been numerically proof for some particular cases |3| . The stability of phase singularities has 
been studied locally [4]. However, definitive conclusions on the behaviour of the singular 
filaments have not been obtained. An alternative approach to this problem could come from 
the study of non-local mechanisms that allow certain topological invariants of the organizing 
centres to be conserved. Preliminary studies on topologically non-trivial field configurations 
have appeared j^, 6|. Examples are the work by Berry and Dennis j3] on phase singularities 
in the Helmholtz equation or studies on the stability of ball lightning j^, ^ . Here we use a 
new approach to the stability of phase singularities in excitable media. We define the helicity 
of the excitable medium and we study its meaning in relation to the total entanglement of 
the system. We consider the time behaviour of the FitzHugh-Nagumo model, showing that 
the persistence of entanglement depends strongly on the conservation of the helicity. This 
observation could be of some utility to develop new methods of controlling the appearance 
of rotors by acting on the system boundaries in experimental situations. 

Excitable media describe in good approximation some properties of certain chemical reac- 
10], cardiac arrhythmias juj, etc. The simplest mathematical models for propagation 



tions 



in three-dimensional excitable media include two state variables u and v that satisfy the 
equations 

— = -f{u,v) + D^V^u, 
at e 

dv 

— = eg{u,v) + D,Vh. (1) 

Here u is the excitation variable, v is the inhibition variable and e is a small parameter. At 
one particular point in the spatial domain T> in which the variables u and v are defined, both 



quantities evolve in time according to (in general, nonlinear) functions f{u,v) and g{u,v) 
respectively. The coupling between close points in the medium occurs due to diffusion terms 
with coefficients and D^. When the system parameters have specific range of values, the 
excitation of a point of the system propagates as a shockwave, with a propagation velocity 
given by the value of e and the diffusion coefficients 12 1. 

In order to study the topology of singular filaments in excitable media, we define a 
vector field in which the field lines coincide with the intersections of level surfaces of u with 
level surfaces of v, since these intersections include the phase singularities that organize the 
complete medium around them. The linking number is a measure of the extent to which 
the field lines of a divergence-free vector field curl themselves around one another, i.e. of 
the helicity of the field as was defined by Moffatt pi5j] in 1969. 

Consider the situation of an excitable media, given by a pair of real scalar fields {u, v) 
defined in a three-dimensional spatial domain V, that satisfy the system of equations (Q. 
Our aim is to describe the unusual stability of certain configurations from a new global 
point of view, taking into account the linkage of the curves obtained from the intersections 
of level surfaces of u with level surfaces of v. Suppose that we are interested in a physical 
situation in which u takes values between Umin and Umax, and v takes values between Vmin 
and Vmax- We define new variables U and V from u and v through linear scaling, in such a 
way that U and V satisfy < f/^ + < 1. Note that the level surfaces of u and v coincide 
with the level surfaces of U and V respectively, since the change is linear. Now we define 
the variables p = \/TP~+V^, q = ArcTan{V/U). A complex scalar field that describes the 
excitable medium is then given by 



^ e'" . (2) 



P 

The level curves of (j) are the intersections of the surfaces of constant u and the surfaces of 
constant v at any time. We now define a vector field given by 

= ^iii^. (3) 

27ii (1 + 00) 

The field lines of Q coincide with the level curves of (p by definition. In equation ([3]), z is 
the imaginary unit and (f) is the complex conjugate of (p. The particular definition ([3]) has 
an interesting bonus. If the scalar is a map (p : S"^, these maps yield knots and can 

be classified in homotopy classes characterized by the integer value of the Hopf index H[(j)), 



which gives a measure of the hnkage or entanglement of the level curves of the map. Since 
the vector field fl defined by equation ([3]) is divergence- free, a vector potential ^ can be 
found so that = V x The Hopf index H{(j)), that is a topological invariant of the map 
(f) : —>■ S"^ , can be then written as the integral 



This integral is known in fluid and plasma physics to be the helicity of the vector field fi, a 
global measure of the linkage of the force lines of fl. Consequently, if we define the helicity 
of the configuration, at any time, as in equation (jlj), this quantity will inform us about the 
global linkage of the intersections of the level surfaces of u and v at any time and, more 
important, if this quantity is conserved during the evolution of the system given by equations 
([T]), then the g'lobal topology of the initial configuration is preserved, constituting a strong 
source of stability of the system. In general, will not correspond to a real map from to 
S"^, so its helicity will not be equal to a Hopf index. This may happen if, for example, the 
domain V is not infinite but a box with finite edges. Then the value of the helicity will be a 
real number instead an integer one, but its meaning is always related to the global linkage 
of the intersections of the level surfaces of u and v. 

An interesting observation can be noted here on the vector field defined from the {u, v) 
configuration through equations (El [3]): it is parallel to the vector field Vu x Vv, whose 
maximum value is used by many authors to detect the vortex that organizes the medium 
{gI. In terms of the global scheme described in this work, the explanation of this fact is 
related with high values of the helicity density in regions where the density of linked lines 
is also high. 

Suppose that a particular initial configuration has been given, in which the initial helicity 
has a non-zero value. When the system evolves in time according to equations ([H), due to 
the presence of diffusion there will be reconnections of the lines given by intersections of 
level surfaces, and the value of the helicity will change in general. However, as we will 
see, there are situations in which the helicity remains constant, or it changes very slowly 
with time compared to the characteristic time of the system. In these situations, the non- 
zero value of the helicity reflexes an unusual stability of the conflguration that can explain 
important numerical, or even experimental, observations. The time variation of the helicity. 
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^ ' ' * ■ Uiv X — (5) 



from equation (jl]), is 

dt Jg " \^"" dt J 
Here, S is the boundary of the three-dimensional domain V and ujv is a unit vector orthog- 
onal to the surface S at each point. The integral ([5]) has to be computed on the boundary 
of the domain. In equation ([5]) we have 



dt 27ri (1 + 00) V^^ 9t 

Expression ([5]) can be used to investigate in what cases the helicity is constant during the 
evolution of the system and how it changes in other cases. A common feature is that the 
conservation of the helicity depends almost only on the boundary conditions, so that the 
stability of the system can be perturbed numerically by acting only on the surfaces of the 
medium. This observation could be of some utility in experimental situations. 

Now let us examine some typical boundary conditions. First, if the domain is the complete 
space and the fields are taken initially so that they are zero at infinity, then the vector 
field d'^/dt will always be zero at the boundaries and the helicity will be conserved in 
time according to equation ([5]). This is the case of the electromagnetic knots in vacuum 



HQ, 
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19| . A similar situation happens if Dirichlet boundary conditions are 



imposed in all the boundaries of a finite box provided the scalars u and v are smooth functions 
of space and the time evolution is also smooth. However, in numerical computations of 
excitable media, Neumann (null flux) and periodic boundary conditions are mostly used. 
As an example, consider that the domain P is a grid in in which the spatial coordinates 
(x, y, z) are confined to the range —L < x,y, z < L, L being a certain length, and suppose 
that Neumann boundary conditions are applied to the x and y directions, and periodic 
boundary conditions are applied to the z direction. In this example, there will not be loss 
of helicity through the z direction, but helicity will be lost through the x and y directions 
in an amount that depends on the parameters of the model that one is computing. Writing 
d"^ /dt in terms of u and v through equations ([2l[6]), the term in equation ([5]) is proportional 
to uat X [dtv'Vu — dtuWv). Taking into account the contributions of both the faces z = —L 
and z = L in which the field is periodic, this is equal to zero, so that the helicity is conserved 
in the directions in which periodic boundary conditions are applied. In the faces in which 
Neumann boundary conditions are applied, the term x {dtvVu — dtuS/v) is not zero but 
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it depends on the parameters of the model because dtu and dtv depend on them through 
the evolution equations ([1]). 

We will give an example of the helicity conservation in the FitzHugh-Nagumo model. 
This model set a paradigma which allows a geometrical explanation of important biological 
phenomena related to neuronal excitability and spike-generating mechanism. It have been 
extensively studied and it was conjectured that persistent solutions called organizing 
centres might exist in three dimensions for excitable media in which two-dimensional vortices 
are embedded into three dimensional space forming knotted vortex rings. In order to prove 
the existence of these solutions, a theoretical framework based on local analysis involving 
effective models of short-range repulsive forces between vortex cores was proposed, but only 
certain limiting cases of slight curvature and twist of vortex lines had partial results 
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20|. Finally to address the fundamental issue of the existence of stable knots, numerical 



investigations were performed. Here we will analyse again the FitzHugh-Nagumo model but 
at the light of the global stability provided by helicity conservation. The FitzHugh-Nagumo 
equations are given by 

= - M - mV3 -v)+ V\ , 
at e ^ ' 

^ = e iu + /3-^v) . (7) 

Here u represents the electric potential and v the recovery variable associated with membrane 
channel conductivity. We choose the constants appearing in equation (JTj) to have the values 
6] £ = 0.3, P = 0.7 and 7 = 0.5. We discretized the equations ([7]) with finite differences on a 
cubic domain of size 2L with a uniform cubic grid of spacing h. For the Laplacian operator 
we use a second order accurate finite difference approximation which is symmetrical up to 
third order. We have to keep in mind the stability criteria about t < h^/2D when setting 
the spatial mesh. 

In order to test the conservation of helicity we enforce u = —1.03279, v = —0.66558 at 
the boundaries, which are the equilibrium values of the system. We take initial conditions 
with nonzero helicity. 



2xz + yjx'^ + y'^ + z^-l) 
u{x,y,z,0) = Ai — 2 I 1^2 "^'^ ' 



( n\ \ - x{x^ + y^ + z'^-l) 
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Here, Ai = a/2 and A2 = l/-\/2 respectively, which are used to cover the range of the 
excitation- recovery loop in {u, v) space for the ordinary differential equation part of 
equations ([7j). Now for testing the conservation of helicity, we will solve the system and 
plot the intersection of level surfaces. Each level surface of constant u can intersect another 
surface of constant f in a curve, few curves or an empty set. If helicity is conserved and 
the initial conditions are chosen in such a way that intersection of level surfaces are knotted 
with entanglement or linking number different from zero, we can find as time evolves that 
there must be intersection curves linked all the time. 

In Figure [1] we plot the evolution of two level surfaces for a domain size L = 5 and 
we have used 100 grid points in each direction. They correspond to the values u = —0.7 
and V = —0.1. Those surfaces like them have a non-trivial intersection. One can use a 
marching cubes algorithm [2l| to get the intersections. In Figure [2] we have plotted the 
intersection of few pairs of {u,v) level surfaces at different times. We can see that, as we 
have shown theoretically, they are linked. It is possible to find linked curves in all the instant 
of times as far we have enough numerical precision. The system eventually will decay to the 
equilibrium values fixed by the boundary conditions, but it will do that keeping the linking 
number constant. In conclusion, we have presented a new global approach to investigate how 
organizing centres in excitable media persist or disappear in three dimensions. In particular, 
we have defined the helicity of an excitable medium as a quantity that describes the global 
linkage of the intersections of the level surfaces of the scalar fields in the medium. We have 
studied how the helicity is conserved or lost through the walls of the medium and shown 
that these behaviours are dominated by the boundary conditions, so the distortion of these 
conditions could lead to the disappearance of the structures. 

We thank S. Betelii for discussions on parts of this research. The authors thank support 
from the Spanish Ministerio de Educacion y Ciencia under project ESP2007-66542-C04-03. 
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& (3 

FIG. 1: Evolution of surface levels. The figure shows the evolution of two level surfaces, u = 
-0.7 (red) and v = -0.1 (green). From left to right and top to bottom, the snapshots correspond to 
0.2, 0.4, 0.6 and 0.8 instants of simulation time. 




FIG. 2: Conservation of entanglement. The figure shows the entanglement due to helicity 
conservation. Each curve results from the intersection of u and v level surfaces. It displays few 
curves at times equal 0.3, 0.5 and 0.8 from left to right. The values of the level surfaces are different 
but at any instant of time there is the same number of linked curves. 
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